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ABSTRACT 

We present a new framework to explain the link between cooling and fragmentation in 
gravitationally unstable protostellar discs. This framework consists of a simple model 
for the formation of spiral arms, as well as a criterion, based on the Hill radius, to 
determine if a spiral arm will fragment. This detailed model of fragmentation is based 
on the results of numerical simulations of marginally stable protostellar discs, including 
those found in the literature, as well as our new suite of 3-D radiation hydrodynamics 
simulations of an irradiated, optically-thick protostellar disc surrounding an A star. 
Our set of simulations probes the transition to fragmentation through a scaling of the 
p hysical op a city. This model allows us to directly calculate the critical cooling time 
of lCammid (l200l . with results that are consistent with those found from numerical 
experiment. We demonstrate how this model can be used to predict fragmentation in 
irradiated protostellar discs. These numerical simulations, as well as the model that 
they motivate, provide strong support for the hypothesis that gravitational instability 
is responsible for creating systems with giant planets on wide orbits. 

Key words: hydrodynamics - radiative transfer - methods: numerical ~ planetary 
systems: protoplanetary discs - planetary systems: formation 
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1 INTRODUCTION 

The fragmentation of protostellar discs through gravita- 
tional instability (GI) is a possible mechanism for the for- 
mation of gas-giant planets and brown dwarfs. For a disc to 
be prone to fragmentation, there are generally thought to 
be two criteria that must be satisfied. 

The first criterion is that a disc must be gravitation- 
ally unstab le. This can b e characterized by the Toomre Q 
parameter l|Toomrel Il964l ). which is the result of a linear 
stability analysis for a differentially rotating thin disc; 

(1) 
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In the above, Cs is the sound speed of the gas, Ke is the 
epicyclic frequency (ke = ^, the rotation rate, for Keple- 
rian rotation), G is the gravitational constant, and E is the 
surface density. For low values of Q ~ 1, such as would be 
found in a massive, cold disc, a disc will be gravitationally 
unstable. 

The second criterion for fragmentation is that a disc, in 
addition to being gravit ationally unstable, must also cool 
quickly. iGammi using 2-D shearing-box simula- 

tions, examined the stability of a local patch of a protostellar 
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disc with a simplified cooling prescription: 

t,oaM = P, (2) 

where fcooi is the cooling time and /3 is a constant. The 
cooling time here is defined through the evolution of the 
specific internal gas energy, u: 

du _ u , , 

By performing simulations with difi^erent values for /?, 
the author found the cooling criterion for fragmentation to 
be 



(4) 



Subsequent work, using global simulations, has shown that 
the critical cooli ng time, Pait, depends on the adiabatic in- 
dex of the gas (|Rice. Lodato fc Armitagel [20051. To date, 
the value of the critical cooling time has only been found 

thr ough num e rical e xperiment. 

ICammid (120011 ). as well as iRice et al] (|2005l ). outlined 
a physical argument for the existence of a critical cooling 
time. If GI can be wel l- charac terized by an a- viscosit y mode l 
llShakura fc SunvaevI l|l973h : see iLodato fc Ric3 (|2004 
l2005h for the applicability of this), then a steady-state can 
exist if the viscous heating by GI is balanced by the pre- 
scribed cooling. If this balance is achieved, then the required 



2 Patrick D. Rogers and James Wadsley 



viscosity of the disc is determined by 

1 
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If GI has a maximum a that it can attain, then it also has 
a maximum heating rate. If the prescribed /3 cools the disc 
faster than this maximum heating rate, then no balance be- 
tween heating and cooling can be achieved, and the disc 
fragments. 

The /3-prescription of cooling, however, is a simplified 
model; more generally, /3 will evolve with a in a given disc. In 
a realistic disc, heating and cooling are linked to the physical 
state of the disc. In this work, we consider fragmentation in 
this more reali stic case, wherea s the /3-prescription of cooling 
does not allow iGammi i (l200lD to have done so. 

Futhermore, the applica bility of a sing l e valu e of /3crit 
(or Qmax) is in some doubt. iMeru fc BatTl|201ll ) gave an 
overview of some of the inconsistencies in the literature re- 
garding the use of a single value for the critical cooling time. 
From a set of their own simulations, the authors found that 
a single value for the critical cooling time was not applica- 
ble; instead, they found an empirical relation in which /3crit 
is a function of the distance from the central star, the lo- 
cal surface density, and the stellar mass. Ho wever, it is not 
clear how well this empirical relation holds, as lMeru fc Bate! 
l|2010l) have demonstrated that 3-D simulations investigating 
the critical cooling time have not yet converged numerically. 

In this work, we present a set of 3-D radiation hydrody- 
namic simulations of a massive, optically-thick, protostellar 
disc, unstable near 100 AU, around an A star. Rather than 
using a /3-prescription for the cooling, these simulations in- 
clude realistic heating and cooling of the disc, including cool- 
ing from the disc photosphere and irradiation from the cen- 
tral star. We do, however, vary the cooling rate in this set of 
simulations by scaling the dust opacity table by a constant 
factor. By reducing the opacity (which reduces the cooling 
time for an optically thick disc) over this set, we observe a 
transition from discs that are stable against fragmentation 
to discs that do f ragment; t his is consistent with the cooling 
criterion work of iGammid (|200ll ). 

W e have used results from ICossins. Lodato fc Clark3 
l|2009l) . and from this set of simulations to develop a sim- 
ple, yet detailed, physical model for the fragmentation of 
a gravitationally unstable protostellar disc. In this model, 
spiral arms develop in an unstable dis c on a characteristi c 
scale related to the disc scale height (jCossins et al.|[2009^ . 
The heating rate of the disc from GI is proportional to the 
square of the amplitude of th e surface density variations in 
the disc (|Cossins et al.l [20091 ): as spiral arms become more 
condensed, the heating rate is increased. The cooling rate 
of the disc from photospheric cooling is inversely propor- 
tional to the square of the surface density; as spiral arms 
become more condensed, the cooling rate in the arms de- 
creases. There is therefore a natural scale for the thickness 
of a spiral arm in a gravitationally unstable disc. This scale 
is set by a balance between heating from spiral waves and ra- 
diative cooling. It is worth noting that for faster cooling rates 
(shorter cooling times), this thickness will be decreased. 

A second scale of interest in this model is the Hill radius, 
which, for an object of mass M, is 

1/3 

^^H.ll= ^ . (6) 
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The importance of the Hill radius can be understood within 
the context of planet formation in a disc of planetesimals 
around a star. If a protoplanet embryo has formed in this 
disc, then it of interest to determine the radius over which it 
may further accrete planetesimals. The Hill radius sets this 
embryo's sphere of influence: material within the Hill radius 
is bound to the embryo and will be accreted. In essence, 
material within the Hill radius of an object is dominated by 
that object's gravity, while material outside of the Hill radius 
is dominated by the central star's gravity, which is equivalent 
to the role of the local shear in the Toomre criterion. For 
the purpose of this discussion, we define the Hill thickness 
as twice the Hill radius. 

In this frame work, we can extend the cooling criterion 
of lGammiel l|200ll ) with the Hill criterion for spiral arms. In 

a gravitationally unstable disc, the natural thickness of the 
spiral arms is set by a balance between heating and cooling. 
Fragmentation occurs in this disc if there is a section of 
arm whose natural thickness is smaller than that section's 
Hill thickness. Essentially, if a section of a spiral arm lies 
within its own Hill thickness, then shear will be unable to 
prevent the collapse of the arm, and fragmentation can take 
place. 

In a gaseous disc, pressure can prevent fragmentation 
from taking place. By considering the Hill radius, we have 
not addressed the role that pressure plays in determining 
fragmentation and how it may modify the critical thickness 
of spiral arms necessary for fragmentation to take place. The 
correct determination of this scale requires the solution to 
a stability analysis of a spiral arm in a difi^erentially rotat- 
ing system. Since we do not have such a solution, we have 
chosen to consider the Hill thickness. The analysis of our 
simulations; however, does indicate that the Hill thickness 
is the correct scale to examine. The Hill criterion for frag- 
mentation is consequently an empirical criterion. 

Our model is consistent with the cooling criterion: as 
the cooling time decreases, spiral arms become thinner and 
more over-dense, becoming more likely to reside within their 
own Hill thickness, and consequently more likely to frag- 
ment. With the Hill criterion, however, we have developed a 
more detailed, and more complete, physical picture of frag- 
mentation. This picture can be applied to the general case of 
a disc with physical heating and cooling, or the more specific 
case of a disc with /^-prescription cooling. In fact, it offers 
a means to calculate what the critical cooling time is for a 
given region of a given disc. 

The structure of the paper is as follows. In ij2l we 
overview our numerical methods as well as our set of sim- 
ulations of gravitationally unstable, irradiated protostellar 
discs. In section ^ we give a detailed picture of our model 
of protostellar disc fragmentation and the Hill criterion. In 
addition, we demonstrate the model's consistency with the 
simulations of ij2l In ^ we show that the Hill criterion is 
quantitatively consistent with the cooling criterion and dis- 
cuss the predictive qualities of the model. Finally, in ^ we 
give our conclusions. 
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2 NUMERICAL SIMULATIONS OF 
GRAVITATIONALLY UNSTABLE 
IRRADIATED DISCS 

2.1 Numerical methods 

Our simulations were performed witli the TreeSPH code 
Gasoline (Wadsley, Sta del & Quinn 2004), with the addi- 
tion of radiative t ransfer in the flux-limite d diffusion ap- 
proximation [FLD: [Rogers fc WadslevI l|201lf) ]. As described 
by the authors, FLD is able to model the transfer of en- 
ergy only in regions in which SPH particles reside. Because 
of limited resolution, any SPH representation of a proto- 
stellar disc naturally has two edges, representing the upper 
and lower atmospheres. Radiative cooling from the disc at- 
mospheres is modelled by means of a photosphere bound- 
ary condition: the SPH particles on the "edge" of the disc 
(the edge-particles) are found, robust surface areas (the area 
of the photosphere for which an edge-particle is responsi- 
ble) are calculated using a 2-D SPH estimate, and a plane- 
parallel cooling term is added to the radiative energy equa- 
tion for the edge-particles. The radiative hydrodynamics has 
been tested on a n umber of sta n dard problems, including the 
relaxation test of iBolev et all (|2007l i. which is particularly 
suited to protostellar disc simulations. 

The conditions in the outer regions of discs (roughly 
100 AU and beyond) are expected to be favourable to 
gravitational fragment ation, since th e cooling criterion is 
likely satisfied there ([Rafikovl |2007| ). As pointed out by 
iKratter fc Murrav-Clavl ( 2011 ) , the heating of the outer re- 
gions is expected to be dominated by the irradiation of the 
disc's surface by the central star, rather than by viscous 
heating. Since our simulations focus on fragmentation at 
these large radii, it is therefore fundamentally important to 
account for this heating via irradiation. 

The photo s phere boundary condition of 
[Rogers fc WadslevI l(201lf) offers a straightforward means by 
which this can be done. In addition to the cooling term in 
the specific radiation energy equation for each edge-particle, 
we have added a heating term of 



Dt 



= —a (Ti, 
rui 



(7) 



where = ^ -I- v • V is the co-moving derivative, ^ is the 
specific radiation energy, Ai is the surface area of the edge- 
particle, rrii is the particle mass, a is the Stefan-Boltzmann 
constant, and Tirrad is the temperature of the stellar irradi- 
ation^ 

iKratter. Murrav-Clay fc YoudinI (|2010D u sed t he pas- 
sive flared disc model of Chiang fc GoldreichI |l993), along 
with a stellar evolution model, to determine the equilibrium 
temperature distribution for a disc surrounding a 1.35 Mq 
A star, which they found to be: 



T = 40 K 



R 



70 AU 



-3/7 



(8) 



Since we are not able to treat the super-heated dust layer of 
optical depth r < 1 in our simulations, it is appropriate to 
use this equilibrium temperature distribution as the ri„ad in 
the irradiation heating term, equation (O. In addition, we 
implement a floor of Ti„ad = 20 K to take into account the 
background radiation field. 
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Figure 1. The physical quantities of the initial disc profile. The 
midplane number density (in units of 10^^ cm~^) is given by the 
solid, black line; the midplane temperature (in units of 10 K) is 
given by the red, dashed line; and the Q is given by the blue, dot- 
dashed line. The horizontal, dotted line is a reference for (3 = 1. 



2.2 Initial conditions and input parameters 

The initial, axisymmetric model of 5 x 10^ SPH p articles 
was c reated in a manner similar to that of Sh en et al.l 
(|2010l )- see Figure [T] for the disc properties. The sur- 
face density profile has the form E oc in the region 
of 20-70 AU (there is a smooth increase of E from 
10-20 AU). There is a smooth functional form of E cx 
r-iexp{-41og(i?) [0.5 log(i?) - log(i?™)] / log(7?o/J?„)}, 
with = 70 AU and Ro = 160 AU, from 70-160 AU, 
after which there is a steep drop off of E cx r~^^ . There is 
roughly 0.61 Mq within 200 AU. 

This particular surface density distribution is motivated 
by the initial temperature profile, which is given by the equi- 
librium temperature in equation ((S|. The combination of 
temperature and surface density results in a broad region of 
the disc having an initial Q of roughly unity. 

Using the initial surface density and temperature pro- 
files, the vertical structure satisfying hydrostatic equilibrium 
was calculated iteratively taking into account both the grav- 
ity from the central star and the self-gravity of the disc. We 
have chosen such a stable initial condition, with Q ~ 1, 
to ensure that the transition towards the spiral structure 
caused by GI is smooth and does not suffer from transients 
that are the result of an unstable setup. 

The central star is modelled as a 1.35 Mq sink par- 
ticle with a radius of 10 AU. We use a mean molecu- 
lar weight of 2.3 and realistic Ros seland-mean opacities 
(|DAlessio. Calvet fc Hartmannll 19971 ). Although the code is 
capable of using a consistent treatment of the internal en- 
ergy of molecular hydrogen that takes in to account transla - 
tional, rotational, and vibrational modes (|Bolev et al.ll2007l ). 
to simplify the analysis of our simulations, we use a fixed 
adiabatic index of 7 = 7/5. The scale height is resolved by 
at least three smoothing lengths outside of 25 AU, and the 
Jeans length is resolved until shortly after fragmentation 
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takes place. We define one ORP to be the Keplerian period 
at our fiducial radius of 100 AU, roughly 863 years. 



2.3 Simulations 

We present a set of five simulations, each of which uses the 
initial conditions described above, the only difference be- 
ing the opacity used. Simulation (A, B, C, D, E) has an 
opacity table that is scaled by a constant value of (1/10, 
1/3, 1, 3, 10). Thus, Simulation C has an estimated phys- 
ical opacity for solar metallicity, while Simulations A and 
B have reduced opacities, and Simulations D and E have 
increased opacities. Physic al changes in opacity could be 
the r esult of grain growth (|Birnstiel. DuUemond fc Brauerl 
l2010l). grain evolution via th e passage of spiral arms 



i Podolak. Mayer fc Quinnll201lh . or formation in an envi- 
ronment with a non-solar metallicity. Our goal, however, is 
not to reproduce different physical environments, but rather 
to explore the necessary conditions for gravitational frag- 
mentation. In this context, a simple scaling of the opacity 
table is both sufficient and desirable. 

All five simulations evolved in a similar fashion over 
the first 2.5 ORPs. High mode-number spiral structure de- 
veloped slowly from SPH Poisson noise in each of the discs 
over this time until settling down to a transitioning state of 
two or three spiral arms. The transition from axisymmetric 
initial condition to spiral structure is observed to be smooth, 
with no strong transients. 

The final states of the simulations are shown in Figure 
[5] Simulations C, D, and E have been evolved for roughly 8.5 
ORPs without fragmentation having taken place (although 
strong spiral two-arm over-densities may persist), while Sim- 
ulation A has fragmented with two objects forming, and 
Simulation B has fragmented with one object forming. This 
set of simulations, therefore, demonstrates a transition from 
non-fragmentation to fragmentation, as a function of the 
opacity scaling. 

For a patch of an optically thick disk, the cooling time 
is approximately 

1 1 C?K^2 



47- 1<jT4~ ' "^^^ 

where n is the opacity l|Rafikov||2007h . Hence, the cooling 
time is directly proportional to the opacity, and our set of 
simulations offers a means of exploring the fragmentation 
boundary as a fun ction of co oling time in a manner similar to 
the simulations of lCamrni^ ( |2001 ,) . The difference is that our 
simulations use realistic radiative cooling (even though the 
opacities may be scaled), rather than /^-prescription cooling. 

We know from the cooling criterion that reducing the 
cooling time (by reducing the opacity) will eventually lead 
to fragmentation. Thus, that Simulations A and B fragment 
is consistent with this picture. However, it is possible to use 
these simulations to better understand why exactly fragmen- 
tation takes place. 

Figure [3] shows the five simulations at roughly the same 
time, shortly before fragmentation took place in Simulations 
A and B. The difference in the structure of the five discs at 
this time offers evidence of a detailed description for why 
Simulations A and B fragment, while Simulations C, D, and 
E do not. As can be observed, as the cooling time decreases 
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Figure 2. The final states of the simulated discs: surface density 
plots of Simulations E, D, C, B, and A are shown from top to 
bottom. As the opacity scaling is reduced from the physical value, 
fragmentation occurs. The discs are shown at respective times of 
8.5, 8.5, 8.5, 3.5, and 2.9 ORPs. 
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(as the opacity decreases), spiral arms in a disc become tiiin- 
ner and more over-dense and this makes fragmentation more 
hkely to occur. 



3 FRAGMENTATION MODEL AND THE HILL 
CRITERION 

We present a model of spiral arm fragmentation in grav- 
itationally unstable discs, based on the observation from 
the simulations of the previous section that reduced cooling 
times lead to thinner arms, which are more likely to frag- 
ment. This model can be broken into two components: the 
first is a model for the (roughly) steady-state spiral struc- 
ture in an unstable disc; while the second is a criterion for 
the fragmentation of these spirals. Many of the details of our 
model are empirical in nature: we h ave used results from the 
simulations of Cossins et al] ||2009| ). as well as our own set 
of simulations in determining some of the important param- 
eters. 



3.1 Spiral structure 

We begin with a model for the spiral structure that results 
in a gravitationally unstable disc. We consider a patch of 
an initially axisymmetric disc that will develop spiral struc- 
ture, such as our initial condition for the simulations of the 
preceding section. This patch is located at some distance, R, 
away from the central star and is of a radial extent lo, with 
a characteristic surface density Sq. GI acts on the scale of 
Iq to collapse mass radialljQ, resulting in the formation of a 
spiral arm of thickness h and characteristic surface density 

El. This process is demonstrated in Figured 

W hat is the appropriate scale for IqI I Cossins et al.l 
(|2009l ) performed a number of simulations of marginally sta- 
ble (Q ~ 1) discs and found that, from a radial Fourier 
transform of these discs, the dominant radial wavenumber 
was typically 



1 



TvGT, 



(10) 



Therefore, we expect the scale of spiral arm formation to 
be lo — 2nH. We have tested that this is consistent with 
our own simulations. Figure [S] shows the results of a radial 
Fourier transform of Simulation B at a time of 2.5 ORPs. 

How many spi ral arms are likely to form in our disc? 
Numerical studies (jLodato fc Rice|[200i , l2005h have shown 
that as the disc-to-star mass-ratio increases, marginally un- 
stable discs show fewer spiral arms. Our simulations are 
of quite massive discs, with a disc-to-star mass-ratio of 
Md/M* = 0.45. This high disk mass is necessary for the 
disc to have Q ~ 1 near 100 AU for our realistic irradiation 
temperature, and results in a typical arm number of m = 2 
or 3. 

The number of arms in a disc is likely the result of swing 
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Figure 3. The simulated discs before fragmentation: surface den- 
sity plots of Simulations E, D, C, B, and A are shown from top to 
bottom. As the opacity scaling is reduced, the spiral arms become 
thinner and more over-dense. The discs are shown at a time of 
2.5 ORPs. 



^ The simulations of ^show tightly- wound spiral structure with 
a typical winding angle of j ~ 10°; to first-order, a purely radial 
collapse is a fair approximation. 
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Figure 4. The creation of a spiral arm: a local patch of radial 
extent Iq (left) in an axisymmetric disc collapses radially to form 
a section of a tightly- wound spiral arm (right) of radial extent l-i . 
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Figure 5. Radial Fourier analysis of Simulation B at 2.5 ORPs: 
the amplitude of the radial mode is given as a function of the 
wavenumber times the local scale height for a number of radii. 
Each amplitude curve is normalized by the peak amplitude for 
th at radius. The Fouri er analysis is carried out using the method 
of lCossins etal] l l200gt) . 



amplification, with significant amplification of a mode, m, 
requiring the swing amplification parameter, 

< 3 l|Binnev fc Trem aine"2008'). From the 



to satisfy 1 < X, 
above, we can see that X, 
star mass-ratio: 



roughly scales witli tlie disc-to- 



Xm oc M^/Ma. 



(12) 



Tlius, for low disc-to-star mass-ratios, only higli-order modes 
will satisfy 1 < Xm < 3, wfiile for iiigh disc-to-star mass- 
ratios, only low-order modes will. 

What is the steady-state thickness, l^, of the newly 
formed spiral arm? We posit that this scale is the result 
of a balance between heating of the disc through the spiral 
waves, and radiative cooling. Assuming that the spiral den- 
sity wave deposits a fixe d fraction, e, o f its en ergy into the 
disc per dynamical time, ICossins et al.l (|2009l ) showed that 



the heating rate per unit mass from spiral arms can written 
as 



where 



M 



kcs 



and M 



kCs 



(13) 



(14) 



are the radial phase Mach number and the Doppler-shifted 
radial phase Mach number, and fip is the pattern speed. 

As outlined by the authors, the pattern speed, and 
hence the Mach numbers, can be calculated from the dis- 
persion relation for a finite-thickness disc: 

27vG\k\ 



m (ilp — Q) 



1 + \k\H' 



(15) 



if the radial and azimuthal wavenumbers are known. From 
their simulations, the authors found a relatively constant 
value of e « 0.2 (see their Figure 15). 

In an irradiated disc, there is additional heating from 
the stellar irradiation, so that 



0+ = 



^ + 



(16) 



The spiral over-density can be calculated by assuming 
that some fraction, /, of the total mass per length within lo 
is compressed into the spiral arm, of thickness h: 



fh 



(17) 



The heating from equation (|16p is balanced by radia- 
tive cooling, for which the cooling rate per unit mass, using 
equation is 



4crT^ 



(18) 



where we have used it = C^/ljij — 1)]. Setting the above 



cooling rate equal to the heating rate of equation (|16p . and 
using the other information in this section, as well as the 
initial axisymmetric properties of our disc, leaves us with 
an equation with only one unknown: the thickness of the 
spiral arm, h, in our patch of interest (assuming that we 
know the proper temperature for radiative cooling, T, in 
the spiral arm. This will be elaborated upon in Q. 



3.2 Determining fragmentation 

Once we know the steady-state thickness of the spiral arm, 
Zi, we can use equation ([BJ to calculate the Hill radius for 
this section of the arm using: 



i^Hill ~ 



30.'^ 



2 1 1/3 



(19) 



If the section of arm has a thickness satisfying 
li/{2Hiim) < 1, then the section lies within its own Hill 
thickness. In the absence of pressure forces, this means that 
the section is bound, as the tidal force from the central star 
(manifest as rotational shear) is less than the self-gravity of 
the section. Once the section is bound, fragmentation oc- 
curs. Conversely, if h / {2Huiu) > 1, then the section of arm 
is not bound and fragmentation does not occur. 
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Figure 6. The Hill criterion for spiral arm fragmentation: if a 
section of spiral arm lies within its own Hill thickness, then that 
section of arm is free to collapse and fragmentation takes place. If 
a section of spiral arm lies outside of its own Hill thickness, then 
shear stabilizes the arm and fragmentation does not take place. 



The Hill thickness tells one about the ability of shear 
to prevent the fragmentation of the arm and is therefore 
expected to be an important scale. However, in comparing 
the radial thickness of the arm to the Hill thickness we are 
ignoring the role of pressure, despite strong radial pressure 
gradients present across the arm. It is therefore reasonable 
to expect that the critical thickness for arm fragmentation 
may be modified from the Hill thickness. There are no strong 
pressure gradients along the arm (the azimuthal direction) 
though, so it is likely that fragmentation occurs in this di- 
rection. In this case, the Hill thickness may in fact be the 
critical scale. 

Determining the correct scale for fragmentation requires 
a detailed calculation of the stability of a spiral arm ac- 
counting for differential rotation. In the absence of such a 
calculation, we posit that the correct scale to consider is 
indeed the Hill thickness. As described below, the results 
from our simulations are consistent with this. The Hill cri- 
terion for fragmentation, demonstrated in Figure [B] is thus 
empirically based. 

3.3 Consistency with simulations 

An analysis of the spiral arms formed in the simulations of 
ij2] shows that their thicknesses, and stability, are consistent 
with the Hill criterion for fragmentation. In Figure [71 we 
show two examples of this analysis to illustrate this consis- 
tency. Our analysis focuses on the surface density of a radial 
slice of the disc (with a typical angular width of 5°). Over 
this slice, a spiral arm is evident as a large over-density. We 
find that arms are often asymmetric; consequently, we de- 
termine a thickness for an arm by fitting each side of the 
arm (with respect to the radius of highest S, -Rp^ak), with a 
Gaussian of the form 



T.0 + e' 



(20) 



where Eo is the value of the surface density adjacent to the 
arm. The thickness of the arm is taken to be h = 2(6ioft -f- 
feright), and the mass of the section of arm is determined 
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Figure 7. Comparison of a spiral arm's thickness to its Hill thick- 
ness. A radial slice of Simulation A (Simulation E) is shown in 
the top (bottom) panel. The surface density of the radial slice 
is given by the black line, while the azimuthally averaged sur- 
face density is given by the blue, dot-dashed line. As described 
in the text, the spiral arm (the large over-density) is fit (the red, 
dashed line) and the arm thickness is found (the vertical, red 
dashed lines). The mass within this arm thickness is computed 
to determine the Hill thickness (the vertical, blue dashed lines). 
Consistent with the Hill criterion, Simulation A (Simulation E) 
has an arm which falls inside of (outside of) its Hill thickness and 
fragmentation is (is not) observed shortly thereafter. 



using a numerical evaluation of 

/-Rpcak-|-26right 
e{R)R'£{R)dR, (21) 

where Q{R) = h/R is the angular extent of the section of 
arm. 

Figure [7] (top) shows this analysis for an arm in Simu- 
lation A, shortly before fragmentation. Consistent with the 
Hill criterion, the arm thickness is less than the Hill thickness 
and therefore fragmentation is expected to occur in this arm; 
indeed, this arm fragmented a short time after the time-step 
used for this analysis. In contrast. Figure [7] (bottom) shows 
this analysis for an arm in Simulation E, which never frag- 
mented. Consistent with the Hill criterion, the arm thickness 
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is greater than the HiU thickness and therefore fragmenta- 
tion is not expected. 

We use these two examples here to illustrate the con- 
sistency of the fragmentation criterion with the simulations 
performed. More generally, we have found that lower opacity 
discs have arms that are consistently smaller with respect to 
their Hill thickness in comparison with higher opacity discs. 
In each of the cases where fragmentation takes place (this 
occurs in the reduced opacity discs) , the arms that fragment 
are observed to fall within their own Hill thickness shortly 
before fragmentation. Indeed, all arms that are observed to 
lie within their own Hill thickness fragment 

3.4 Comments regarding the model 

The simplicity of the model described above is desirable, as 
it gives a straightforward, physical picture of the formation 
of spiral arms in a gravitationally unstable disc, as well as 
the physical criterion that determines whether or not those 
spiral arms fragment. However, unstable discs do exhibit a 
great deal of complexity; here, we discuss this complexity 
and comment on its implications for our model. 

We have described the steady-state thickness of a single 
spiral arm as being the result of a balance between heating 
and cooling. However, the simulations show that the spiral 
structure in the disc evolves with time: the number of arms 
in the disc is not constant, nor is the over-density of each 
arm. However, an analysis of the spiral arms in these sim- 
ulations shows that the radial sound crossing time of these 
arms is generally much less than the radiative cooling time. 
As a consequence, the arms are able to adjust very quickly 
to perturbations, so that contractions are quasi-static. The 
important point is that the arm thickness responds to the 
balance between heating and cooling. 

The Hill criterion describes the fragmentation of a spi- 
ral arm, but it does not necessarily determine whether or 
not this fragmentation leads to the formation of a long- 
lived object; this also depends on the cooling of the frag- 
ment and the complex environment of the disc in which 
the initial fragmentation takes place. Simulation C, for ex- 
ample, demonstrates an instance of "failed fragmentation," 
as shown in Figure |8] One of the spiral arms appears to 
have fragmented; however, the resulting over-density is only 
short-lived: it shortly thereafter collided with the next spiral 
arm, without surviving. 

The Hill criterion describes the formation of a fragment 
based on the inability of shear to stabilize a section of a 
spiral arm. Further collapse occurs on the radiative cool- 
ing timescale of the fragment. If this timescale is long, then 
the fragment may still be quite diffuse, and easily disrupted 
by collisions with subsequent spiral arms. Indeed, Simula- 
tion B, which generally has shorter cooling times than Sim- 
ulation C because of its opacity scaling, shows a fragment 
which formed, but then subsequently collided twice with spi- 
ral arms; in contrast to Simulation C, this object survived, as 
observed in Figure [2] Fragmentation, therefore, can be well 
characterized by our model; however, whether or not frag- 
ments survive also depends on the complex non-linear inter- 

^ There is one exception: the spiral arm of Simulation C that 
collided with a fragment, described in the next section. 




100 AU 



Figure 8. An instance of "failed fragmentation" in Simulation 
C. In this example, a fragment, the over-density within the black 
circle, is observed; however, it is not long-lived. Shortly after this 
time-step, the fragment collided with the spiral arm and did not 
survive. 



actions between collapsing fragments and the spiral struc- 
ture in the disc. 

As described, the fragment of Simulation C was dis- 
rupted through a collision with the subsequent spiral arm. 
This resulted in a strong compression of the spiral arm; in 
fact, the compression was strong enough that the arm was 
observed to lie within its own Hill thickness. Nevertheless, 
the arm did not fragment. This does not conflict with the 
Hill criterion because in this instance, the arm was not in a 
near steady-state. Since the timescale for the collision was 
much shorter than the sound-crossing time of the arm, the 
arm could not adjust. As a result, the increased over-density 
of the arm lead to an increased heating rate, see equation 
HHI, and a reduction of the cooling time, see equation @. 
As a consequence of the imbalance between heating and 
cooling, pressure forces caused the arm thickness to expand 
on roughly the sound-crossing time, with the result that no 
fragmentation took place. 



4 CONSISTENCY WITH THE CRITICAL 

COOLING CRITERION AND PREDICTIVE 
ABILITY 

The physical model developed in the preceding section can 
be used to examine whether or not fragmentation is likely 
to take place in a disc. In this section, we demonstrate that 
the predictions of this model are consistent with previous 
numerical results of discs evolved using the /3-prescription 
of cooling, as well as the results of our suite of simulations 
discussed in ij2l Specifically, we analyze the initial condition 
of our simulated disc described in [J2]and adopt general val- 
ues of /o = 2-kH, e = 0.2, MM = 1 and / = 1. These values 
will be correct to within 0(1), but will likely have variation 
depending on the physical properties of a disc. 
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Figure 9. The critical cooling time (black line), as calculated 
for our disc initial condition. The horizontal blue, dot-dashed line 
represents the value as previously determined by numerical ex- 
periments, while the red, dashed line shows the Q profile. The 
vertical, red lines show the extent of the unstable region with 
Q < 1.4. As can be seen, the value for /3crit in this region is 
consistent with the value from numerical experiment. 



4.1 Calculating the critical cooling time 

The critical cooling time, /3crit, for a marginally stable disc 
with 7 = 7/5 ha s been found to b e /3crit = 12 from numeri- 
cal experiments [Rice et aP (j2005l ): with the caveat that nu- 
merical convergence has not been clearly demonstrated]. If 
we adopt a heating rate without irradiation (consistent with 
the aforementioned simulations) , given by equation (|13p and 
balance this heating with a /3-prescription cooling rate, given 
by 



7(7-1)/?' 



(22) 



then we can solve for the cooling rate that results in a certain 
arm thickness: 



/3 = 



eMM'y (7 — 1) 



3/2 / ofr,'2 \ 1/2' 



(23) 

The critical cooling time can then be computed from the 
above equation by setting the arm thickness to exactly 
match the Hill thickness, h / (2Hum) ~ 1- 

The critical cooling time computed for our disc initial 
condition is shown in Figure|9] As can be seen, we do not find 
a unique value for the critical cooling rate, but rather a value 
that depends strongly on radius (due to the variation of disc 
properties with radius). Importantly, we observe that our 
model of disc fragmentation predicts critical cooling times 
in the unstable region of the disc (where Q ~ 1) that are 
consistent with the results of numerical experiments. 

This consistency is a useful check on our model. In ad- 
dition, it is noteworthy that Figure [9] demonstrates the first 
calculation of the critical cooling time from a physical model 
of fragmentation. Previous estimates of /3crit have come only 
from numerical experiments. 



Figure 10. The spiral arm thickness (black curves) of an irradi- 
ated disc in units of the Hill thickness, as calculated for our disc 
initial condition. The arm thickness is calculated for the range 
of opacities used in the simulations. From the curve of greatest 
h/(2-ffHiii) to the curve of smallest Zi/(2//Hiii) the opacity scal- 
ings are 10, 3, 1, 1/3, and 1/10 the physical opacity (the solid 
black curve). The Q profile is given by the red, dashed curve, 
while the vertical, red lines show the extent of the unstable re- 
gion with Q < 1.4. 



4.2 Predictive ability of the model for irradiated 
discs 

It is useful to check our model against the results of previous 
work using cooling in the form of a /3-prescription. However, 
it is of particular interest to apply the model to the more 
realistic case of an irradiated disc with radiative cooling. 
Without considering GI, an irradiated disc has a natural 
equilibrium state in which the heating of a particular patch 
of disc from stellar irradiation is balanced by the cooling 
of that patch from the radiative cooling of the disc photo- 
sphere. Here, we consider deviations from this equilibrium 
state due to GI. 

Specifically, there is an additional heating of the disc 
from the spiral arms, given by equation (jlSp . which will re- 
sult in an increase in temperature of ST. We consider spiral 
arms in which this excess heating is balanced by a perturba- 
tive radiative cooling; that is, a cooling rate for the dissipa- 
tion of this excess thermal energy. The perturbative cooling 
time for an irradiated disk, when the temperature in the arm 
is comparable to the irr adiati on temperature (as we find for 
the discs in is (e.g. iKratter et al.ll201ol ) 



7 



32 7 - 1 crT*,. 



-kT, 



(24) 



Balancing the perturbative heating and cooling, we 
solve for the natural arm thickness for our disc initial con- 
dition from !j2]in units of the Hill thickness. This is done by 
solving the following quartic: 



{C}lt 



B 



: 0, (25) 
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where 



B ■- 
and 

C ■■ 



eclMMQ. 



f 



3k7 V^o^o, 
Once li is calculated, we can calculate the ratio 



2-ffHill 



r 

'1 



2/3 



GEoZo 



1/3 



(26) 



(27) 



(28) 



The results of these calculations, for the range of opacity 
scalings used in our simulations, are shown in Figure [Tol 

We expect discs to fragment for /i/(2HHiii) ^ 1, and the 
trend in Figure[lO]is consistent with this picture: as found in 
our simulations, the increased-opacity discs are less likely to 
fragment (have larger Zi/(2_ffHiii)) than the reduced-opacity 
cases. However, even though the increased-opacity discs are 
not expected to fragment, the arm thickness is expected to 
be within a factor of two of the Hill thickness in the region 
of Q ~ 1. This result is consistent with our simulations, as 
is demonstrated for the arm from Simulation E shown in 
Figure [7| (bottom). This arm is stable, in that it does not lie 
within its Hill thickness; however, its thickness only exceeds 
the Hill thickness by a relatively small factor. 

From our calculation, only the lowest opacity case, cor- 
responding to Simulation A, has h / {2Huiu) < 1, and would 
be expected to fragment; in fact, both Simulation A, and 
the second-lowest opacity case, Simulation B, have shown 
fragmentation. This discrepancy is likely simply the result 
of the choice of parameters used in the calculation. We have 
chosen values for a number of parameters in our model {lo, 
f, e, and MM) that are expected to be correct to within 
0(1); however, the exact values will likely have some vari- 
ation. With an improved understanding of the growth of 
spiral structure, and the heating of spiral arms, the model's 
predictive abilities will be improved. 

In the inner regions of irradiated discs, our applica- 
tion of a perturbative cooling is not expected to be ap- 
propriate. In these regions, the heating is expected to be 
dominated by viscous heating , rather than by irradiation 
l|Kratter fc Murrav-Clavll2011^ ■ which may result in a disc 
temperature that is significantly different from the irradi- 
ation temperature. The breakdown of our calculation can 
be observed as the vertical jumps in /i/(2i?Hiii) observed in 
Figure 1101 such as the jump at 80 AU for the highest opac- 
ity disc. Inside of this jump radius, our calculations do not 
give the correct result. It is possible to calculate the proper 
h/{2Huiu) in this region by balancing the total heating and 
cooling from equations (|16|) and (|18|l . rather than the per- 
turbative variants; however, this requires a knowledge of the 
temperature within the spiral arms. Therefore, we leave an 
analysis of these regions for future work. 



5 DISCUSSION AND CONCLUSIONS 

5.1 Implications for planet formation 

Direct-imaging observations have shown the existence of 
gcis-giant planets at large distances from their host A 
star, including HR 8799b, 7 M.,up at a distance of 68 AU 



(|Marois et alj2008h. and Fom alhautb, 3 Mj^p at a distance 
of 119 AU l|Kalas et all 120081 ). It is difficuh to explain the 
existence of gas-giants at such distances from their host star 
in the core accretion scenario, since the surface densities are 
typically too low to form th e necessary rocky cores within 
the lifetime o f the gas disc (|Dodson-Robinson et al ] |2009l : 
iRafikovHioTH ). However, more investigation is warranted 
in order to determine if such planets can be explained in 
the core- accretion scenario. In comparison, fragmentation 
via GI has been shown to be a viable formation mecha- 
nism at large dista nces from the host star from both the- 
oretical arguments l|Rafikovl [20071 : iNero fc BiorkmanI [20091 : 
iKratter et aL[|2010l '). as well as from numerical simulations 
with radiative transfer (|Bolevll2009l : Istamatellos et al.ll201ll : 
lBossll201ll ). 

The particular set of 3D radiation hydrodynamic simu- 
lations presented here were designed to investigate fragmen- 
tation at large radii (~ 100 AU) around A stars. At these 
distances, heating from stellar irradiation is expected to be 
the dominant heating source; we have includ ed irradiation 
using the Ti„ad expected for a 1.35 M0 A star (jKratter et al.l 
l2010l i. 

The results of our simulations show that GI can pro- 
duce gravitationally bound objects at large distances from 
the star, given opacities on the low side of the expected 
rang e. Such op acities could be the result of grain-growth 
( Birnstiel et al..,2 010), gr ain ev olution via the passage of spi- 
ral arms iPodolak etaL[|201ll ). or formation in an environ- 
ment with a no n-solar metallicity (HR 8799 is roughly 1/3 
solar metallicity, iMarois et al.ll2008l ) . Our simulations do not 
take these physical mechanisms into account, but rather use 
a simple scaling of the opacity table. 

Although our simulations do not have the resolution to 
follow the evolution of bound objects as their central den- 
sities run away, it is interesting to consider the objects at 
the end- state of our simulations, shown in Figure [5] At the 
end of Simulation A (with an opacity scaled by 1/10), there 
are two brown dwarfs of masses 21 Mjup and 15 Mjup, at 
respective distances of 62 AU and 95 AU; while at the end 
of Simulation B (with an opacity scaled by 1 /3) , there is one 
brown dwarf of mass 40 Mj^p at a distance of 95 AU. Nei- 
ther the masses nor the distances of these objects represent 
their final state: all of the objects are accreting mass and 
migrating inwards at the end of the simulation. 

We conclude that GI in unstable discs can produce 
brown dwarfs at large distances from A stars. We have, of 
course, only shown fragmentation for a single surface den- 
sity and temperature profile. It is of interest to investigate 
a greater region of the parameter space with numerical sim- 
ulations in order to explore the po ssibility of low - mass com- 
panions suc h as th ose observed bv iMarois^ et al.l (|2008l ) and 
iKalas et al.1 l|2008l ). 

5.2 Physical model of fragmentation 

We have presented a new framework to explain the link be- 
tween cooling and fragmentation in protostellar discs. This 
framework consists of two components. The first is a sim- 
ple model for the formation of spiral arms, in which the 
thickness of a spiral arm is set by a balance between heat- 
ing (through gravitational instability and irradiation) and 
radiative cooling. The second is a criterion for fragmenta- 
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tion: spiral arms that have a natural thickness smaller than 
their Hill thickness fragment, resulting in objects that may 
survive to become gas-giant planets or brown dwarfs. 

This model of fragmentation is based on results from 
ICossins et all (|2009l ') as well as our suite of 3-D radiation 
hydrodynamics simulations of gravitational instability in an 
irradiated, optically-thick protostellar disc surrounding an 
A star. By reducing the opacity scaling, and consequently 
the cooling time, over the set of simulations, we have pro- 
duced a suite that demonstrates the transition from non- 
fragmentation to fragmentation. From an analysis of these 
simulations, we have found that the critical scale for deter- 
mining fragmentation is roughly the Hill thickness: those 
spiral arms that are found to fragment lie within twice their 
Hill radius, while those spiral arms that do not fragment ex- 
tend beyond their Hill thickness. In the future, it would be 
of interest to have a robust calculation of the critical scale 
for fragmentation from a stability analysis of a spiral arm in 
a differentially rotating system. 

In comparison to the critical cooling time picture, our 
model of fragmentation is a more detailed, and more general, 
physical picture of fragmentation that is applicable to discs 
with realistic heating and cooling, in addition to discs with 
/3-prescription cooling. Indeed, by coupling the Hill criterion 
to our simple model of spiral arm formation, heating, and 
cooling using a /3-prescription, we have been able, for the 
first time, to calculate P^^it . We find that there is not a single 
value for /3crit, but that it depends on the local properties of 
the disc; in addition, our calculation is consistent with the 
value determined by numerical experiment. 

We have also demonstrated how this model can be used 
to predict fragmentation in irradiated discs with radiative 
cooling. Applying the model to the initial condition of our 
simulated disc, for the various opacity scalings used, yields 
predictions that are consistent with the results of our sim- 
ulations. An improvement in the predictive abilities of the 
model depends on a better understanding of several param- 
eters that describe the formation and heating of the spiral 
arms. 

This model has been developed in the context of pro- 
tostellar discs; however, it may also be of use in the con- 
text of star formation in a disc near the Galactic centre 
l|Levin fc BeloborodovlboOSl ). as well as star-cluster forma- 
tion in optically-thick starburst galaxies, such as Arp 220. It 
is beyond the scope of this work to consider these systems, 
and we leave these considerations to future work. 

In this work, we have considered isolated discs; that 
is, the effects of accretion from the surrounding envelope 
were ignored. However, accretion is expected to play an im - 
portant role in grav itationally unstable discs [BolevI l|2009l ): 
iKratter et al] l|20ld )]. since it is accretion that will push the 
mass of the disc towards being sufficient for instability to set 
in, keep it unstable despite mass-transport, and contribute 
to heating. In future work, we intend to investigate the ef- 
fects of accretion on the stability of protostellar disks, in the 
context of our model of fragmentation. 
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